# TODO: Add comment
# 
# Author: wuping
###############################################################################

samplingQuick.angleOrder=function(sampleMethod,populationData){
	if(sampleMethod@composite){
		n=sampleMethod@n
		sampleMethod@n=round(n/2)
		pointToevent=eventToevent=sampleMethod
		pointToevent@focusType="point"
        eventToevent@focusType="event"
        d1=single.angleOrder(pointToevent,populationData)
		d2=single.angleOrder(eventToevent,populationData)
		return(list(d1,d2))
	}else{
		return(single.angleOrder(sampleMethod, populationData))
	}
}

#TODO add a overlap situation
single.angleOrder=function(sampleMethod,populationData){
	
	#prepare
	#amplify the sample size for each focus point
	if(!sampleMethod@overlap){
		if(sampleMethod@focusType=="point")
			spt=sampleMethod@n*3
		else{
			spt=sampleMethod@n*3
			nnn=length(populationData@availible)
			if(spt>nnn){
				spt=nnn
			}
		}
	}else{
		spt=sampleMethod@n+10
	}
	
		
	endPoints=cbind(populationData@x,populationData@y)
	kmax=sampleMethod@k*sampleMethod@q*6
	if(kmax>dim(endPoints)[1])
		kmax=dim(endPoints)[1]-1
	d=NA
	redo=TRUE
	delna=numeric()
	reptime=0
	
	while(redo & reptime<=sampleMethod@stopLimit){
		
		#step1 select focus points
		if(sampleMethod@transectLine){
			dline=transectLine(sampleMethod,populationData)
			xrange=dline@xrange
			#transect line only appliable for point to event
			x=runif(spt,xrange[1],xrange[2])
			y=dline@slope*x+dline@b
		}else{
			if(sampleMethod@focusType=="point"){
				x=runif(spt,sampleMethod@edge, populationData@plotdim[1]-sampleMethod@edge)
				y=runif(spt,sampleMethod@edge, populationData@plotdim[2]-sampleMethod@edge)
			}else{
				select=sample(populationData@availible,spt,replace=sampleMethod@overlap)
				x=populationData@x[select]
				y=populationData@y[select]
			}
		}
		focusPoints=cbind(x,y)
		
		#step2 find kmax neighborhoods
		kn=get.knnx(endPoints, focusPoints, k=kmax,algorithm="kd_tree")
		#truncate the index and distance matrix
		if(sampleMethod@focusType=="event"){
			kn$nn.index=kn$nn.index[,-1]
			kn$nn.dist=kn$nn.dist[,-1]
		}
		
		#step3 find which sector each selected neighbor is
		if(sampleMethod@q>1){
			
			angle=angleToFocus(endPoints,focusPoints,kn,sampleMethod@q)
			
			sector=t(apply(angle,1,nearestMatch,q=sampleMethod@q,k=sampleMethod@k))
			delna=unique(which(is.na(sector))%%spt)
			delna[which(delna==0)]=spt
			if((length(delna)+1)>=dim(sector)[1]){
				reptime=reptime+1
				next()
			}
				
			if(length(delna)>0){
				sector=sector[-delna,]
				kn$nn.index=kn$nn.index[-delna,]
				kn$nn.dist=kn$nn.dist[-delna,]
				angle=angle[-delna,]
			}
			localspt=dim(sector)[1]    
			#transform to position in a vector
			position=as.vector((sector-1)*localspt+1:localspt)
			#position=position[!is.na(position)]
			#browser()
		}else{
			position=(1:spt)+(sampleMethod@k-1)*spt
		}
		
		#step4 check is there enough non-replicated sample lines
		nleft=length(position)
		nrowleft=nleft/sampleMethod@q
	
		#browser()
		index=(kn$nn.index)[position]
		#uniqueindex=unique(index)
		
		#firstap=match(uniqueindex,index)
		#firstap=firstap%%nrowleft
		#firstap[which(firstap==0)]=nrowleft
		#nuniquerow=length(uniqueindex)
        if(!sampleMethod@overlap){
			uniqueindex=uniquerows(index,sampleMethod@q)
			#print(c(length(position),length(unique(position))))
			uniquerowindex=attr(uniqueindex,"rown")
		}else{
			uniqueindex=index
			uniquerowindex=1:nrowleft
		}
        dim(position)=c(nrowleft,sampleMethod@q)
		uniqueposition=position[uniquerowindex,]
		dim(uniqueindex)=c(length(uniquerowindex),sampleMethod@q)
        nuniquerow=length(uniquerowindex)
		
		if(nuniquerow<sampleMethod@n){
			redo=TRUE
			reptime=reptime+1
		}else{
			redo=FALSE
			#position.unique=position[match(unique(index),index)]
			#get the distances 
			#选择时并不考虑这些距离是否需要凑足一个中心点的几个q，另一中取法是去取中心点。
			#select=sample(position.unique,sampleMethod@n*sampleMethod@q)
	        
			#dim(uniqueposition)=c(nuniquerow,sampleMethod@q)
			
			
			#position=position[firstap,]
			#nrowleft=length(firstap)
			selectrow=sample(1:nuniquerow,sampleMethod@n)
			if(sampleMethod@q==1)
				select=as.vector(uniqueposition[selectrow])
			else
			    select=as.vector(uniqueposition[selectrow,])
			d=(kn$nn.dist)[select]
			
			dim(d)=c(sampleMethod@n,sampleMethod@q)
			
			#browser()
			###################
			#visual testing code
	        if(sampleMethod@visual){
				plot(endPoints,xlab="x",ylab="y")
				if(length(delna)>0)
				focusPoints=focusPoints[-delna,]
				#browser()
				focusPoints=focusPoints[uniquerowindex,]
				
				focusPoints=focusPoints[selectrow,]
				if(is.null(dim(focusPoints)))
					dim(focusPoints)=c(1,2)
				points(focusPoints,pch=4,col=2,lwd=2)
				index=as.vector(uniqueindex[selectrow,])
				#print(c(dim(angle),dim(kn$nn.index)))
				#print((angle)[select])
				#print(c(length(index),length(unique(index))))
				
				toPoints=endPoints[index,]
			    points(toPoints,pch=19,col=3)
				
				fx=rep(focusPoints[,1],sampleMethod@q)
				fy=rep(focusPoints[,2],sampleMethod@q)
				arrows(x0=fx,y0=fy,x1=toPoints[,1],y1=toPoints[,2],length=0.1)
				for(i in 1:dim(focusPoints)[1]){
					fp=focusPoints[i,]
					b=tan(seq(0,2*pi,2*pi/sampleMethod@q))
					a=fp[2]-b*fp[1]
					for(j in 1:sampleMethod@q)
					    abline(a=a[j],b=b[j],lty=2,col="grey")
				
				}
				#browser()
			}
			
		}
	}
	if(reptime == sampleMethod@stopLimit)
		return(FALSE)
	else
	    return(d)
}


uniquerows=function(index,q){
	do=TRUE
	rown=1:(length(index)/q)
	while(do){
		nrowleft=length(index)/q
		tb=table(index)
		dbn=as.numeric(names(tb))[tb>1]
		if(any(tb>1)){
			firstap=match(dbn,index)
			firstap=firstap%%nrowleft
			firstap[which(firstap==0)]=nrowleft
			un=unique(firstap)
			dim(index)=c(nrowleft,q)
			index=as.vector(index[-un,])
			rown=rown[-un]
		}
		do=length(index)!=length(unique(index))
	}
	attr(index,"rown")=rown
	return(index)
}

nearestMatch=function(x,q,k){
	x0=1:q
	for(i in 1:k){
		re=match(x0,x)
		if(i < k){
			x[re]=NA
		}
	}
	return(re)
}

angleToFocus=function(endPoints,focusPoints,kn,q){
	dim.nn=dim(kn$nn.index)
	#dim.nn[2]=dim.nn[2]-1
	nn.index=as.vector(kn$nn.index)
	point=endPoints[nn.index,]
	dist=as.vector(kn$nn.dist)
	ac=(point[,2]-focusPoints[,2])/dist
	ac[ac>1]=1
	ac[ac<(-1)]=-1
	#range from -pi/2 to pi/2
	angle=asin(ac)
	#divide into two groups
	negative=angle<0
	#find points might in 2 or 3 image
	cng=point[,1]<focusPoints[,1]
	
	im2=!negative & cng
	im3=negative & cng
	angle[im2]=pi-angle[im2]
	angle[im3]=abs(angle[im3])+pi
	angle[angle<0]=angle[angle<0]+2*pi
	
	angle=try(as.numeric(cut(angle,breaks=seq(0,2*pi,2*pi/q),include.lowest=T)))
		
	dim(angle)=dim.nn
	return(angle)
}


